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I 1^ Abstract. We investigate the dynamics of vortices in repulsive Bose-Einstein 

fi -s , condensates in the presence of an optical lattice (OL) and a parabolic magnetic 

•-H ' trap. The dynamics is sensitive to the phase of the OL potential relative to the 

■ • magnetic trap, and depends less on the OL strength. For the cosinusoidal OL 

potential, a local minimum is generated at the trap's center, creating a stable 
equilibrium for the vortex, while in the case of the sinusoidal potential, the vortex 
J^ , is expelled from the center, demonstrating spiral motion. Cases where the vortex 

is created far from the trap's center are also studied, revealing slow outward- 
spiraling drift. Numerical results are explained in an analytical form by means 
of a variational approximation. Finally, motivated by a discrete model (which 
is tantamount to the case of the strong OL lattice), we present a novel type of 
vortex consisting of two pairs of anti-phase solitons. 
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The experimental realization and theoretical studies of Bose-Einstein condensates 
(BECs) pP have led to an explosion of interest in the field of atomic matter waves and 
their nonlinear excitations, including dark fl\ and bright [3] solitons. More recently, 
two-dimensional (2D) excitations, such as vortices 0] and vortex lattices |3, were 
considered and realized experimentally. Other nonlinear states, such as e.g., Faraday 
waves jQ, ring dark solitons and vortex necklaces 0, stable solitons and localized 
vortices in attractive BECs trapped in a 2D optical lattice (OL) [H], and even stable 
solitons supported by an OL in a repulsive BEG jH] were also predicted. 

Vortices, in particular, are worth studying not only due to their significance as 
a fundamental type of coherent nonlinear excitations, but also because they play a 
dominant role in the breakdown of superflow in Bose fluids |10[lll| . The theoretical 
description of vortices in BECs can be carried out in a much more efficient way than 
in liquid He |12| due to the weakness of interactions in the former case (which, in 
addition, is tunable |13p. These advantages explain a large volume of work regarding 
the behavior of vortices in BECs, which has been recently summarized in a review 

The subject of the present paper is the dynamics of vortices under the action 
of the OL; to the best of our knowledge, this problem is considered in this work for 
the first time. An OL potential is generated as an interference pattern by laser beams 
illuminating the condensate, in particular in ID and 2D cases [151 1161 [T7lll8[ll9| . In 
the 2D case, the OL potential assumes the form (in dimensionless units) 

VoL{x,y) - Vo [cos'^ikx + cf,) + cos^ {ky + (j))] , (1) 

where Vq is the strength of the OL, which is measured in units of the recoil energy 
Er — h? /2m\^,^^^ (i.e., the kinetic energy gained by an atom when it absorbs a 
photon from the OL), where Aiascr is the laser wavelength, h is Planck's constant, 
m is the atomic mass, k is the wavenumber of the OL, and is a phase-detuning 
factor [the obvious possibility to remove (j) by means of the diagonal shift, {x, y) — > 
(x — (j^l^iV ~ 4>/k)i is ruled out by the presence of the magnetic-trap potential, see 
below] . The wavenumber k of the OL can be experimentally controlled by varying the 
angle between the counter-propagating lasers producing the interference pattern j2()j . 
As is well known, the effective 2D GP equation applies to situations when the 
condensate has a nearly planar ("pancake") shape, see, e.g., [2] and references 
therein. Accordingly, vortex states considered below are not subject to 3D instabilities 
(corrugation of the vortex axis (14|) as the transverse dimension is effectively 
suppressed. We base the analysis on the Gross-Pitaevskii (GP) equation, written 
in harmonic-oscillator units |22| . 

iut = —Au + \u\ u + V{x, y) u, (2) 

where u = u{x,y,t) is the 2D wave function, and the external potential is 

V{x,y)^^n^ix^+y^) + VoL{x,y), (3) 

which includes the isotropic magnetic trap 1 and OL potential ^. Notice that 
the dimensionless parameter Q = cur/ujz in Eq. JSjl, where ujr,z are the confinement 
frequencies in the radial and axial directions respectively, is assumed to be 51 ^ 1. In 
most studies of vortices, angular momentum imparted by stirring of the condensate is 
assumed (see, e.g., [331 1^ and references therein), as typically this is the way in which 
vortices are generated in experimental settings ^ . Here we assume that a vortex has 
been generated, but then the stirring ceases. Alternatively, the vortex may be created 
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by means of the phase-engineering technique j25| . The positive sign in front of the 
nonhnear term in Eq. ^ imphes that we consider a repulsive condensate. 

We first present numerical findings for the case of a "regular" vortex; subsequently, 
the observed dynamics is explained by means of a variational approximation. We then 
proceed to examine a new type of a vortex, which is a robust bound state of two pairs 
of TT-out-of-phase pulses, as suggested by results known for 2D dynamical lattices. In 
most cases, we fix the strength of the magnetic-trap potential to be fl^ ~ 0.002. This 
value is relevant to a ®^Rb condensate of radius 25/xm, containing 4.6 x 10^ atoms in 
a highly anisotropic trap with Ur = 2Tr x 7.5 Hz and ojz = 2tt x 115 Hz. According 
to these values of the physical parameters, in the following results, which will be 
presented in normalized time and space units, the corresponding units are 1.38 ms 
and 1/im respectively. Additionally, as far as the OL parameters are concerned, the 
typical values Vq = 0.5 and k = 1 used in most cases, correspond to O.bEr and to an 
OL wavelength of 6.3/im 

In simulations, the initial vortex configuration is taken, in polar coordinates p 
and 0, as 

u — p(r) exp {iO) utf, (4) 

where p{r) = r'^{Q.M + 0.07r2)/(l + QAlr"^ + 0.07r'') is the radial Fade interpolation 
for the vortex solution to the GP equation without external potential, and wtf = 
^max {0, p — (51^/2) [x^ + j/-^)} is the Thomas-Fermi (TF) wave function for the 
magnetic trap ^, with a chemical potential p. The vortex is placed at the center 
of the magnetic trap (unless otherwise indicated). Simulations were performed by 
means of a finite-difference discretization in a box of the size 60 x 60 (60/xmx60/xm), 
using 200 X 200 points (i.e., the spatial stepsizes are dx — dy — 0.3). Time integration 
was performed by means of the fourth-order Runge-Kutta scheme with the time step 
dt = 0.0025 {dt — 3A5ps). It should be remarked that in the initial stages of the time 
evolution, the initial condition of Eq. Q) "adjusts" itself to the presence of the OL, by 
shedding small amplitude radiation waves. These are dissipated by the an absorbing 
region close to the boundary implemented in the numerical simulation. 

The first case that was examined is the one with </) = in Eq. (^. In this case, 
the vortex is extremely robust, staying at the center during a few hundred time units 
of Eq. (O, which correspond to a few hundred of milliseconds in physical units. As 
can be seen in Fig. ^ where the vortex and its motion are shown for Vq — 0.5 and 
fc = 1, the displacement of the vortex center remains, indefinitely long, as small as 
~ lO^'^. The center of the vortex was located by fitting, close to the local minimum, 
the density jwp to the expression |up = Ax'^ + By^ + Cx + Dy + E, which yields the 
position of the center at {-C/2A, ~B/2D)). 

In Fig. 121 we display the evolution which starts from the same initial configuration, 
but with (j) = 7r/2 in Eq. (^). A distinctly different behavior is observed in this case, 
which is shown for longer times, up to t = 250 (346 ms), to verify that it is a real 
feature. In particular, it is clearly observed in this case that the vortex spirals out 
from the center of the magnetic trap. Note that fine details of the motion indicate 
meandering motion, similar to that exhibited by spiral waves in excitable media j26j . 
even though in the present case the meandering occurs on a smaller spatial scale, 
hence much finer resolution would be needed to clarify it. However, the outward- 
spiraling motion of the vortex center is clearly seen in Fig. 13 Considerably longer 
runs [for times on the order of a few thousand time units of Eq. (j21) , which correspond 
to seconds in physical units] show that the vortex spirals all the way to the edges of 
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Figure 1. The vortex on a BEC of diameter ^ 50^m in tlie cosinusoidal [</> = 
in Eq. Q] OL potential with Vq = 0-5 and fc = 1. The top panel shows the 
contour plot of the density at t = 100 (138 ms). The bottom left panel is a cut 
of the same density profile along x = 0, while the bottom right one shows the 
motion of the vortex center for < t < 100, the initial position being marked by 
a star. Notice the scale (10~^, or 1 nm in physical units) of very weak jiggling of 
the vortex, which thus stays practically immobile at the origin. 






Figure 2. The same as in Fig.Q but now the potential is sinusoidal, = tt/2. 
To clearly demonstrate the dynamics, the top and bottom left panels show the 
density profiles at a larger time than in Fig. Q t = 250 (346 ms). The bottom 
right panel depicts the motion of the vortex center for < t < 250 [positions of 
the vortex center at t = 100 (138 ms) and t = 200 (276 ms), respectively, are 
indicated by the star and circle]. 



the TF cloud where it eventuaUy decays (see also comments below). 

In Fig. O we examine the case with — tt/A in Eq. ^. In this case also, 



the 
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Figure 3. The case of Vq = 0.5 and <j) = '''/4 in Eq. 0. The vortex moves 
rapidly away from the center. The snapshots show the density profile (top panel) 
and its horizontal and vertical cuts (along x = 1.5 and y = —1.5, respectively) of 
the field at t = 100 (138 ms). 



vortex moves outward. In fact, it is observed that it moves away from the trap center 
much faster, as at i = 100 (138 ms) it is already located near {x,y) = (1.5, —1.5) (cf. 
Fig. [21), corresponding to a couple of /xm from the trap center. 

Qualitatively similar conclusions, but at different time scales, have been obtained 
for different values of the OL strength Vq. For instance, in the potential with cj) = 7r/2, 
Vo — 1.5 (and k = 1), the vortex spirals out again (similar to the case with Vq = 0.5 
displayed above), the time necessary for it to reach the distance of 0.1 from the origin 
being three times as large as in the case of Vb = 0.5 (details are not shown here). 
This can be explained by the fact that, if Vq is larger, the vortex has to move through 
a "rougher" energy landscape, hence it becomes more difficult for it to "find its way 
out". 

These findings can be qualitatively understood in terms of an effective potential 
which governs the motion of the vortex; in particular, the potential has a minimum 
and maximum at the center of the magnetic trap in the cases oi (j) = and 4> = 7r/2, 
respectively. To present this explanation in a mathematically tractable form, we resort 
to a variational approximation. As in Ref. j27| . we use the following ansatz for the 
vortex, inspired by its linear analog namely the first excited state of the quantum 
harmonic oscillator: 

u(a;,y,t) = i?(t)ro(t)exp[-r2(t)/6(t)] e«^«(*), (5) 

where rg(t) ee (x - xo(t)f + (y - 2/o(t)) ^, and LpQ{t) = i&n-'^[{y - yo{t))/{x ~ Xf){t))]. 
Note that the ansatz carries a vortex-like structure centered at (a;o,j/o)j and it 
qualitatively emulates the initial waveform Q which was adopted in the numerical 
simulations. The stronger the nonlinearity, the less accurate the approximation offered 
by Eq. |SJ| is; however, as the ansatz (O was found to always provide a qualitatively 
correct description of the phenomenology, we use it to represent the vortex. 

Similarly to the calculations reported in Refs. |211 ES|, one can easily deduce 
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(additionally using the norm conservation) that the evolution of B{t) and b{t) is, to 
the leading order, negligible, therefore we set them to constant values, B « B{t) 
and b « b{t). Then, the substitution of the ansatz jSJ) into the Lagrangian of the 
GP equation Q with the potential l|3Il leads, up to constant factors, to the effective 
Lagrangian 

^ = 2 (*o + Vo) - K.ff (xo, yn), (6) 

where the net effect of the OL and magnetic trap is combined into an effective 
potential, 

T4ff (x, y) = Q{^) [cos{2kx) + cos{2ky)] + ^n^ {x^ + y^) , (7) 

where Q{(j)) is given by a rather cumbersome expression; we will actually need the 
values 

g(0) = -^ {bk^ - 2) exp(-6fcV2), 
8 



Q (I) = -g(0), and 



0. 

Equations © and |7l suggest that the motion of the coordinates xq and j/o of 
the vortex core resembles the motion of two uncoupled oscillators. Actually, this is a 
straightforward generalization of a result that can be obtained in the corresponding 
ID case: in that case, the counterpart of the vortex is a dark soliton, whose equation 
of motion in the presence of the potential l|2Jl can be derived by means of the adiabatic 
perturbation theory for dark solitons in BECs; see e.g., [SOj and references therein. 
This approach, in the presence of the one-dimensional magnetic trap and OL, yields 
an effective potential 



vn^)=\vo 



\"' z\k' 

6 



cos(2fca;) + -fl'^x'^. 



In fact, Eq. Q is a generalization of this expression. The effective potential from Eq. 
(Q is depicted in Fig.^as a function of x and y for the cases ^ = and (j) = tt/2, where 
we set b — 1 [this value was chosen, comparing the size of the numerically obtained 
vortex with that implied by the ansatz ©J. 

In our experiments fc = 1 was used, and thus, in the case </> = 0, a vortex 
positioned at the center of the magnetic trap is at a local minimum of the potential 
(O, therefore the vortex stays in this position. If slightly perturbed around this stable 
fixed point, the vortex will perform small amplitude oscillations. On the contrary, in 
the case of </> = tt/2, the same point {xq, yo) — (0, 0) is unstable (as a saddle), and as 
a result the vortex moves away from it. These conclusions are in full accord with the 
numerical experiments, see Figs. ^a-ndEl 

Normally, in a Hamiltonian mechanical system the motion from an unstable saddle 
point occurs along the corresponding unstable separatrix. However, it is well known 
pn 1221 133 that a spatially-periodic setting may give rise to resonant effects, which 
induce an effective dissipation, due to emission of radiation waves from a nonlinear- 
wave state (vortex, in the present case). This effective dissipation, as is known 
P^ in^ 03 ESI J causes the vortex to spiral outward (making the saddle point look 
like an unstable spiral), which is what we indeed observe in the simulations, see Fig. 
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Figure 4. The effective potential Veff{x,y), obtained in the variational 
approximation [cf. Eq. Q] , responsible for the motion of the vortex subject to the 
2D optical lattice and magnetic trap. Top panel: the cosinusoidal optical lattice 
{(j> = 0) produces a minimum in the effective potential at the origin. Bottom: the 
sinusoidal optical potential {(p = tt/2) gives rise to a maximum of the potential 
at the origin. The right panels display more clearly the combined effect of the 
optical lattice and magnetic trap (here we use D.^ = 0.005, instead of fl^ = 0.002, 
which was used in the computations). The parameters are 6 = 1, A; = 1, Vb = 0.5. 



[3 Note that in other cases where the Hamihonian is not positive definite (such as 
e.g., Korteweg-de Vries - type models), radiative losses may destabilize an equilibrium 
position which is otherwise stable jSZ]- Finally, in the case of (/) = 7r/4, the center of 
the vortex is not originally at an equilibrium position of the potential of Eq. Q , hence 
its drift is manifested faster, as is indeed seen in Fig. |3| 

Now we turn to a different case, in which the vortex was created far from the 
center of the magnetic trap (the results are displayed only for (f) — 0, as it was found 
that the value of 4> does not significantly affect the results in this case). The motion 
of the vortex is displayed in Fig. where its outward-spiraling, due to the effective 
dissipation induced by the optical lattice, is obvious. This spiral motion is combined 
with jiggling induced by the potential energy landscape. Thus, we conclude that a 
vortex seeded in the periphery of the BEC cloud will slowly drift towards the edge of 
the cloud, where it will eventually decay into other excitations |35| . 

Finally, we present a novel type of vortex, which is quite different from the one 
considered above. Motivated by recent investigations of vortices on 2D discrete lattices 
PHI EH) I we initialized a real field configuration which is shown in the top panel of Fig. 
El The configuration consists of two up-down (dipolar) pulse pairs, each featuring a 
phase shift of tt, hence a total phase shift along a contour surrounding these pulses is 
27r, which corresponds to unit vorticity. In discrete lattices, numerically exact stable 
stationary solutions of this type exist [SH]- In spatially uniform continuum models, 
such stationary solutions cannot exist, but Fig. shows that they may be sustained 
by the OL. After shedding some radiation, the configuration reaches a nearly steady 
state, although residual oscillations are observed. Notice, once again, that waves 
radiated away due to the oscillations are absorbed by the dissipative boundaries 
of the integration domain (which were imposed to avoid artificial reflection). This 
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Figure 5. The top and middle panels show the vortex configuration at t = and 
t = 750 (~ 1 s). The bottom panel shows the trajectory of the vortex from t = 
up to t = 750 (the two positions are marked by the star and circle, respectively). 



dynamically stable, lattice-motivated configuration represents, to the best of our 
knowledge, a novel vortex-like structure which is particular to the system equipped 
with OL: without the lattice potential, the two 7r-out-of-phase sets of the in-phase 
pulses would separate due the pulse-pulse interactions (which are repulsive in the case 
of opposite "parity" pulses and attractive in the case of same parity ones; see e.g., 
Uni)- However, the effective potential exerts a "local force" on the pulses, making it 
possible to trap them together in a stable configuration of this type. 

In conclusion, we have studied basic properties of vortices in the 2D Gross- 
Pitaevskii equation with the optical lattice (OL) and magnetic trap. The results 
crucially depend on the phase of the OL relative to the parabolic magnetic trap. 
Depending on the phase, it is possible to trap the vortex at the center of the trap, 
or, on the contrary, to expel it, which is readily explained in terms of an effective 
potential for the vortex derived by means of the variational approximation. In the 
latter case, the vortex moves along an unwinding spiral, which is explained by the 
effective dissipation due to emission of small-amplitude waves. We also considered the 
case in which the vortex was initially created far from the center of the trap, for which 
a precessing motion with a slow outward-spiraling drift was observed and qualitatively 
explained. Finally, a new species of a stable vortex, specific to the system equipped 
with the optical lattice, was found in the form of a bound state of two 7r-out-of- 
phase soliton pairs. The existence of this vortex is also explained by the OL-induced 
potential. 

Naturally, it would be of interest to further examine structures which are inspired 
by features particular to periodically modulated models, and to observe how they 
appear and disappear with variation of the OL strength. Such studies are currently 
in progress. 
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Figure 6. The top panel shows the contour plot of the real part of the initial 
configuration consisting of two pairs of pulses which are n out of phase (hence 
the total phase change, as one goes around the configuration is 2tt). The pluses 
indicate (centers of) the "up-pulses" , whereas the minus signs the "down-pulses" 
of the initial configuration. The middle and bottom panel show the contour plot 
of the square modulus of the wave function at t = 100 and t = 200 respectively. In 
the latter, we have verified that the total phase shift around the pulse quartet is 
27r, hence the vortex retains its topological charge. In this case, Vb = 2, fc = 0.5, 
(i = 0, and C = 0. 
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